SUBROUTINE standarderrors(np,p)
  USE commonvars
  IMPLICIT NONE
  INCLUDE 'mpif.h'

  INTEGER, INTENT(IN)  :: np
  REAL(8), INTENT(IN)  :: p(np)
  INTEGER              :: i, j, N_points, info, ii, tt, Npoly, np_p
  INTEGER, PARAMETER   :: nop = 10, Nfinest = Nfines-1
  REAL(8),ALLOCATABLE  :: dsim(:,:,:), fval_vec(:), hh(:), mag(:), work(:)
  REAL(8)              :: deps, da, h, newp(np+nop), dg(np,Nmom), dm(nop,Nmom), dgdp(np+nop,Nmom), fval, fval_init, fval_min, pp(np+nop), sumbetadf, uu, stderr(np)
  REAL(8)              :: Vtemp(np,Nmom), VV(np,np), Varcov(np,np), dsim_orig(Nmom), Vpf(4,4), Vw(6,6), Vop(nop,nop), inv_dg(Nmom,np), Vtemp_op(Nmom,nop), &
                          VV_op(Nmom,Nmom), VarCov_temp(Np,Nmom), VarCov_op(np,np), TotVarCov(np,np), VarCovPol(np,np), Id_temp(np,np)
  REAL(8)              :: diffmom(Nmom,1), temp(Nmom,Nmomdata), temp1(4,4), temp2(6,6), pfine_temp, pfine_temp1, &
                          cdf_fine(Nfinest), true_cdf_fine(Nfinest), frac_fine_temp, frac_finet(Nfinest), &
                          qcons_w_riv, prob_a_temp, prob_a_temp1, cdf_a(Nab), true_cdf_a(Nab), a_discr_temp, beta, res_prod_func
  INTEGER              :: lwork, iparam, imom, ipiv(np), count_statsig
  CHARACTER            :: uplo ='U'

  EXTERNAL utilfun, emaxfun, simulatemun, DGELSD, DGETRF, DGETRI
  
  flag_lowerfval = 0

  N_points = 4
  ALLOCATE(dsim(Np+nop,N_points,Nmom),fval_vec(N_points),hh(N_points),mag(np+nop))
  
  !Machine precision
  DEPS = 1.0d0
  
2 DEPS = DEPS/2.0d0
  DA = DEPS + 1.0d0
  
  IF(DA .GT. 1.0d0) THEN
     GO TO 2
  ELSE
     DEPS = ABS(2.0d0 * DEPS)
  END IF

  IF (myrank == 0) write(*,*) 'DEPS', DEPS

  ii = 0

  ii = ii + 1
  ! theta is the relative value of public consumption
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        theta = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        theta = factor_param(ii)*p(ii)
     END IF
  ELSE
     theta = 0.5d0
  END IF
  ! Used in the computation of standard errors
  pp(ii) = theta
  
  ii = ii + 1
  !theta_v is the relative value of public consumption for voters; we switched to a model with only one type in terms of relative value
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        delta_fine = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        delta_fine = factor_param(ii)*p(ii)
     END IF
  ELSE
     delta_fine = 0.1d0
  END IF
  ! Used in the computation of standard errors
  pp(ii) = delta_fine
  
  ii = ii + 1
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        cost_run_frac = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        cost_run_frac = factor_param(ii)*p(ii)
     END IF
  ELSE
     cost_run_frac = 0.1d0
  END IF    
  pp(ii) = cost_run_frac
  
  ii = ii + 1
  ! upower
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        upower = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        upower = factor_param(ii)*p(ii)
     END IF
  ELSE
     upower = 1d0
  END IF
  pp(ii) = upower
  
  ii = ii + 1
  ! Electoral appeal parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        eta = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        eta = factor_param(ii)*p(ii)
     END IF
  ELSE
     eta = 1d0
  END IF
  pp(ii) = eta

  ii = ii + 1
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        sigma_run = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        sigma_run = factor_param(ii)*p(ii)
     END IF
  ELSE
     sigma_run = 1d0
  END IF
  pp(ii) = sigma_run

  ii = ii + 1
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        sigma_steal = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        sigma_steal = factor_param(ii)*p(ii)
     END IF
  ELSE
     sigma_steal = 0d0 ! 0.05142432d0
  END IF
  pp(ii) = sigma_steal

  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_ie = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_ie = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_ie = 0d0
  END IF
  pp(ii) = alpha_ie

  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_nst = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_nst = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_nst = 0d0
  END IF
  pp(ii) = alpha_nst
  
  ii = ii + 1
  ! Electoral appeal parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_ea = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_ea = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_ea = 0d0
  END IF
  pp(ii) = alpha_ea
  
  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_q = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_q = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_q = 0d0
  END IF
  pp(ii) = alpha_q
  
  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        alpha_st = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        alpha_st = factor_param(ii)*p(ii)
     END IF
  ELSE
     alpha_st = 0d0
  END IF

  pappeal(1) = 1d0
  
  ii = ii + 1
  ! Effect of being audited and having stolen on election probabilities parameter
  IF (ii <= nparam) THEN
     IF (flag_sa == 0) THEN
        zeta = factor_param(ii)*(ubp(ii)*EXP(p(ii))+lbp(ii))/(1d0+EXP(p(ii)))
     ELSE
        zeta = factor_param(ii)*p(ii)
     END IF
  ELSE
     zeta = 0d0
  END IF
  
  !We use a log-normal distribution for the fine
  !mean of the fine distribution
  mu_fine = 0.0942693d0 + delta_fine
  !standard deviation of the fine distribution
  sigma_fine = 0.2837051d0

  !The following calculations produce the discrete approximation for the probability of drawing one of the discrite values, frac_fine(:), and one of the corresponding discrete values, pfine(:)
  !Given the parameters of the fine distribution we derive the support points that gives the best approximation (part B in Kennan)
  ! WE USE Nfinest = Nfines-1 to compute the distribution of fines; we then add one more point that corresponds to the case of no conviction with corresponding probability
  DO i = 1, Nfinest
     ! cdf values at the fixed support points
     true_cdf_fine(i) = (2d0*DBLE(i)-1)/(2d0*DBLE(Nfinest))
     !We compute the inverse of the uniform at true_cdf_fine(i); we consider a uniform defined on [0,mu_fine+sigma_fine*E[stealing in data]]
     CALL normal_01_cdf_inv(true_cdf_fine(i),frac_fine_temp)
     frac_finet(i) = EXP(sigma_fine*frac_fine_temp + mu_fine)
  END DO

  !We compute the probabilities at the chosen points that best approximate the log-normal distribution (part A in Kennan)
  CALL normal_01_cdf((LOG(frac_finet(1)) - mu_fine)/sigma_fine,pfine_temp)
  DO i = 2, Nfinest
     CALL normal_01_cdf((LOG(frac_finet(i)) - mu_fine)/sigma_fine,pfine_temp1)
     cdf_fine(i-1) = (pfine_temp + pfine_temp1)/2d0
     pfine_temp = pfine_temp1
  END DO
  cdf_fine(Nfinest) = 1d0

  ! This corresponds to the case of no conviction
  frac_fine(1) = -1d0
  frac_fine(2:Nfines) = frac_finet(:)

  !We compute the probability of the discrete fine values in our grid
  ! pfine(1) corresponds to the probability of no conviction
  pfine(1) = 1d0 - prob_conv
  pfine(2) = cdf_fine(1)
  DO i = 2, Nfinest
     pfine(i+1) = cdf_fine(i) - cdf_fine(i-1)
  END DO
  pfine(2:Nfines) = pfine(2:Nfines)*prob_conv
  
  IF (flag_bootstrap == 0) THEN
     IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf_fine', true_cdf_fine(:)
     IF(myrank == 0) write(*,'(a)') ''
     IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf_fine', cdf_fine(:)
     IF(myrank == 0) write(*,'(a,4f20.10)') 'frac_fine', frac_fine(:)
     IF(myrank == 0) write(*,'(a,5f20.10)') 'pfine', pfine(:), SUM(pfine(:))    
  END IF

  !The ability of mayors is drawn from LOG-N(mu_a,sigma_a)
  !Mean of the ability distribution, we normalize the constant in the production function to 0 and assign the entire constant in the IV estimation to mu_a
  res_prod_func = 3.0858715d0  !3.078571d0
  mu_a = -3.516575d0 + res_prod_func ! -0.2941902d0 !-0.32856629d0 !(ubp(1)*EXP(p(1))+lbp(1))/(1d0+EXP(p(1))) !-.53632455d0
  !standard deviation of the ability distribution
  sigma_a = 0.4391229d0  !0.39916454d0 !0.44185212d0 ! .35204222d0 !.3438509d0 !.4730158d0 ! .48566405d0
  !Mean of the ability distribution in the second term = mu_a + beta' - beta, where beta' and beta are the estimated constants in the production functions for the first and second term
  beta = mu_a !-.2941902d0 !-0.32856629d0
  
  !parameters of the production function for public consumption                                         
  !we use the same constant in both terms
  !power parameter for public inputs
  alphaq(1) = 0.1941907d0  !.1526144d0 !0.1518577d0
  !power parameter for private inputs
  alphaq(2) = 0.4339804d0  !.429166d0 !0.3725852d0
  !linear term in education
  alphaq(3) = 0d0  !0.0430096d0 !0.0340758d0
  !We normalize the constant to be equal to 0
  alphaq(4) = beta - mu_a
  !flexible polynomial in population
  alphaq(5) = 0d0 !6.66d-06 
  alphaq(6) = 0d0 !-2.26d-10
  alphaq(7) = 0d0 !1.46d-15
  alphaq(8) = 0d0 !-3.75d-21
  alphaq(9) = 0d0 !3.16d-27

  DO i = 1, 2
     pp(nparam + i) = alphaq(i)
  END DO
!!$  DO i = 4, Nparq-1
!!$     pp(nparam + i) = alphaq(i+1) 
!!$  END DO
  pp(nparam + 3) = mu_a - res_prod_func
  pp(nparam + 4) = sigma_a
    
  !Polynomial of order n in population
  Npoly = Nparq - 4
  polypop = 0d0
!!$  DO imun = 1, Nmun
!!$     DO i = 1, Npoly
!!$        polypop(imun) = polypop(imun) + alphaq(4+i)*popsize(imun)**DBLE(i)
!!$        !if (myrank ==0) write(*,'(a,2i4,4f20.10)') 'polypop', imun, i, polypop(imun), alphaq(4+i), popsize(imun), popsize(imun)**DBLE(i)
!!$     END DO
!!$  END DO
  
  !The following calculations produce the discrete approximation for the probability of drawing one of the discrite values, prob_a(:), and one of the corresponding discrete values, a_discr(:)
  !Given the parameters of the ability distribution we derive the support points that gives the best approximation (part B in Kennan)
  DO i = 1, Nab
     ! cdf values at the fixed support points
     true_cdf_a(i) = (2d0*DBLE(i)-1)/(2d0*DBLE(Nab))
     !We compute the inverse of the standard normal at true_cdf_a(i)
     CALL normal_01_cdf_inv(true_cdf_a(i), a_discr_temp)
     a_discr(i) = sigma_a*a_discr_temp + mu_a
  END DO
  
  !We compute the probabilities at the chosen points that best approximate the log-normal distribution (part A in Kennan)
  CALL normal_01_cdf((a_discr(1) - mu_a)/sigma_a,prob_a_temp)
  DO i = 2, Nab
     CALL normal_01_cdf((a_discr(i) - mu_a)/sigma_a,prob_a_temp1)
     cdf_a(i-1) = (prob_a_temp + prob_a_temp1)/2d0
     prob_a_temp = prob_a_temp1
  END DO
  cdf_a(Nab) = 1d0

  !We compute the probability of the discrete fine values in our grid
  pability(1) = cdf_a(1)
  DO i = 2, Nab
     pability(i) = cdf_a(i) - cdf_a(i-1)
  END DO

  !We take the exp of a_discr(i) to make it a log-normal
  a_discr(:) = EXP(a_discr(:))

  IF (flag_bootstrap == 0) THEN
     IF(myrank == 0) write(*,'(a)') ''
     IF(myrank == 0) write(*,'(a,3f20.10)') 'true cdf ability', true_cdf_a(:)
     IF(myrank == 0) write(*,'(a,3f20.10)') 'discrete ability', a_discr(:)
     IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf ability', cdf_a(:)
     IF(myrank == 0) write(*,'(a,3f20.10)') 'prob discrete ability', pability(:)
  END IF
  
  !parameters of the wage function for past mayors
  !constant
  alphaw(1) = -3.086317d0
  !education
  alphaw(2) = 4.096486d0
  !age
  alphaw(3) = 1.035222d0
  !age^2
  alphaw(4) = -.0975276d0
  !audit*ln(popsize)
  alphaw(5) = 0d0
  !audit*d_steal*ln(popsize)
  alphaw(6) = 0d0
  !dummy for medium size munic
  alphaw(7) = 1.144047d0
  !dummy for large munic
  alphaw(8) = 3.215137d0
  !standard deviation
  sigmaw = 6.335244d0 !6.464272d0

  DO i = 1, 4
     pp(nparam + 4 + i) = alphaw(i)
  END DO
  DO i = 1, 2
     pp(nparam + 8 + i) = alphaw(6 + i)
  END DO
  
  !Standard deviation of the private input shocks
  sigma_zzpr = 0d0 !(ubp(20)*EXP(p(20))+lbp(20))/(1d0+EXP(p(20)))

  !Standard deviation of the shocks to the decision to run
!!$  sigma_run = 1d0 ! 1d0 !0.1d0 ! (ubp(19)*EXP(p(19))+lbp(19))/(1d0+EXP(p(19)))
  
  ! WE WILL ASSUME LOG PREFERNCES
  !preference parameters                                                                                
  delta = 0d0 ! 2d0 !0d0 ! p(6) if we use log utility for public consumption, otherwise we have to estimate the parameter 
  !gamma = 1 if we use log utility for private consumption, otherwise we have to estimate the parameter 
  gamma = 2d0

  ! WE SHOULD GET RID OF THE MEASURMENT ERRORS
  !Mean of the measurement error in stealing; a positive mean means that it is more likely that a positive clearical mistake (leading to corruption) is made than negative (leading to negative corruption)
  !As a percentage of funds (see simulatemun)
  mu_steal = 0d0 !(ubp(22)*EXP(p(22))+lbp(22))/(1d0+EXP(p(22)))

  !We transform the level of public consumption to take into account the degree of rivalry
  !This is also the variable we output for public consumption
  !If we divide by (popsizeg)**eta, the utility at subsistence level goes to -infinity 
  qcons_w_riv = subqcons_vec(2) + small_no

  subutil = 0d0
  CALL utilfun(subpcons,qcons_w_riv,theta,uu)
  IF (flag_bootstrap == 0 .AND. myrank==0) WRITE(*,*) 'uu,subpcons,subqcons,theta', uu,subpcons,qcons_w_riv,theta
  sumbetadf = 0d0
  DO tt = 1, Tend
     sumbetadf = sumbetadf + betadf**(tt-1)
     subutil(Tend - tt + 1) = uu*sumbetadf
  END DO

  IF (flag_bootstrap == 0 .AND. myrank==0) WRITE(*,*) 'subutil', subutil(:)

  IF (flag_bootstrap == 0 .AND. myrank == 0) THEN
     write(*,*) ''
     write(*,*) 'param', p(:)
     write(*,*) 'electoral probability: alpha_ea, alpha_st, alph_q, alpha_ie', alpha_ea, alpha_st, alpha_q, alpha_ie, alpha_nst
     write(*,*) 'pappeal', pappeal
     write(*,*) 'mu_fine, sigma_fine, delta_fine', mu_fine, sigma_fine, delta_fine
     write(*,*) 'plability', pability
     write(*,*) 'prod fn: pu*ab, pr*ab, educ, const', alphaq(1:4)
     write(*,*) 'prod fn: poly pop', alphaq(5:Nparq)
     write(*,*) 'theta', theta
     write(*,*) 'upower', upower
     write(*,*) 'wages past mayors', alphaw
     write(*,*) 'delta, gamma', delta, gamma
     write(*,*) 'sigma_run, sigma_zzpr', sigma_run, sigma_zzpr
     write(*,*) 'sigma_steal, mu_steal', sigma_steal, mu_steal
     write(*,*) 'cost_run_frac, upower, upower-cost_run_frac', cost_run_frac, upower, upower-cost_run_frac
     write(*,*) 'mu_a, sigma_a', mu_a, sigma_a
     write(*,*) 'beta', beta
     write(*,*) 'eta', eta
     write(*,'(a,11f18.14)') 'param inputs', theta, delta_fine, cost_run_frac, upower, eta, sigma_run, sigma_steal, alpha_ie, alpha_nst, alpha_ea, alpha_q
  END IF
    
  CALL emaxfun()
  
  CALL simulatemun()

  !Store the results of each size of perturbation j for each parameter i!
  dsim_orig(:) = simmom(1:Nmom,1)
  
  IF (myrank == 0) write(*,'(a,11f20.10)') 'dsim', dsim_orig(:)

  CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)
  
  !----------- COMPUTE THE INITIAL FVAL ------------!

  ! Compute the difference in moments
  diffmom(1:Nmom,1) = simmom(1:Nmom,1) - datamom(1:Nmom,1)

  DO i=1,Nmom
     temp(i,1)=SUM(diffmom(:,1)*INV_V(:,i))
     if (myrank == 0) write(*,'(A,1I3,x,A,1F20.8,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6,x,A,1F12.6)') 'i',i,'dfval',diffmom(i,1)*INV_V(i,i)*diffmom(i,1),'simmom',simmom(i,1),'datamom',datamom(i,1),'diff',diffmom(i,1),'V',V(i,i),'weight',INV_V(i,i)
  END DO
  if (myrank == 0) write(*,*) ''

  fval=SUM(temp(:,1)*diffmom(:,1))

  IF (myrank == 0) THEN
     write(*,*) 'initial fval for stderrs, IT HAS NO MEANING BECAUSE THE MOMENTS FOR THE PROD FN PARAMETERS ARE DIFFERENT FROM THE DATA', fval
  END IF

  fval_init = fval
  fval_min = fval
  
  dgdp = 0d0 

  !Array for the 5 point derivative method!
!!$  hh = (/ -2d0, -1d0, 1d0, 2d0 /)
  hh = (/ 1d0, 2d0, 3d0, 4d0 /)

  !Parameter adjusting size of perturbation
  DO i = 1, nparam+nop
     mag(i) = 1d0
  END DO

  mag(6) = 4d0
  mag(7) = 3d0
  
  IF(myrank == 0) write(*,'(A,21F10.3)') 'mag', mag
  IF(myrank==0) write(*,'(a,21F20.10)') 'Params', pp(:)
  
  np_p = np+nop
  
  DO i = 1, np_p
     DO j = 1, N_points

        h = hh(j)*deps**(1d0/3d0)*ABS(pp(i))*mag(i)
        IF(myrank == 0) write(*,'(A,1F12.8,1I3,3F16.8)') 'param,i,hh(j),h,mag',pp(i),i,hh(j),h,mag(i)

        newp(:) = pp(:)
        newp(i) = pp(i) + h
        IF(myrank == 0) write(*,'(A,4F20.12)') 'newparam',newp(i),pp(i),h,pp(i)+h

        theta = newp(1)
        delta_fine = newp(2)
        mu_fine = 0.0942693d0 + delta_fine
        cost_run_frac = newp(3)
        upower = newp(4)
        eta = newp(5)
        sigma_run = newp(6)
        sigma_steal = newp(7)
        alpha_ie = newp(8)
        alpha_nst = newp(9)
        alpha_ea = newp(10)
        alpha_q = newp(11)
        
        mu_a = newp(np+3) + res_prod_func
        sigma_a = newp(np+4)

        !power parameter for public inputs
        alphaq(1) = newp(np+1)
        !power parameter for private inputs
        alphaq(2) = newp(np+2)
        !We normalize the constant to be equal to 0
        alphaq(4) = beta - mu_a

        alphaw(1) = newp(np+5)
        !education
        alphaw(2) = newp(np+6)
        !age
        alphaw(3) = newp(np+7)
        !age^2
        alphaw(4) = newp(np+8)
        !dummy for medium size munic
        alphaw(7) = newp(np+9)
        !dummy for large munic
        alphaw(8) = newp(np+10)

        !The following calculations produce the discrete approximation for the probability of drawing one of the discrite values, frac_fine(:), and one of the corresponding discrete values, pfine(:)
        !Given the parameters of the fine distribution we derive the support points that gives the best approximation (part B in Kennan)
        ! WE USE Nfinest = Nfines-1 to compute the distribution of fines; we then add one more point that corresponds to the case of no conviction with corresponding probability
        DO ii = 1, Nfinest
           ! cdf values at the fixed support points
           true_cdf_fine(ii) = (2d0*DBLE(ii)-1)/(2d0*DBLE(Nfinest))
           !We compute the inverse of the uniform at true_cdf_fine(ii); we consider a uniform defined on [0,mu_fine+sigma_fine*E[stealing in data]]
           CALL normal_01_cdf_inv(true_cdf_fine(ii),frac_fine_temp)
           frac_finet(ii) = EXP(sigma_fine*frac_fine_temp + mu_fine)
        END DO
        
        !We compute the probabilities at the chosen points that best approximate the log-normal distribution (part A in Kennan)
        CALL normal_01_cdf((LOG(frac_finet(1)) - mu_fine)/sigma_fine,pfine_temp)
        DO ii = 2, Nfinest
           CALL normal_01_cdf((LOG(frac_finet(ii)) - mu_fine)/sigma_fine,pfine_temp1)
           cdf_fine(ii-1) = (pfine_temp + pfine_temp1)/2d0
           pfine_temp = pfine_temp1
        END DO
        cdf_fine(Nfinest) = 1d0
        
        ! This corresponds to the case of no conviction
        frac_fine(1) = -1d0
        frac_fine(2:Nfines) = frac_finet(:)
        
        !We compute the probability of the discrete fine values in our grid
        ! pfine(1) corresponds to the probability of no conviction
        pfine(1) = 1d0 - prob_conv
        pfine(2) = cdf_fine(1)
        DO ii = 2, Nfinest
           pfine(ii+1) = cdf_fine(ii) - cdf_fine(ii-1)
        END DO
        pfine(2:Nfines) = pfine(2:Nfines)*prob_conv
        
        IF (flag_bootstrap == 0) THEN
           IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf_fine', true_cdf_fine(:)
           IF(myrank == 0) write(*,'(a)') ''
           IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf_fine', cdf_fine(:)
           IF(myrank == 0) write(*,'(a,4f20.10)') 'frac_fine', frac_fine(:)
           IF(myrank == 0) write(*,'(a,5f20.10)') 'pfine', pfine(:), SUM(pfine(:))    
        END IF

        !The following calculations produce the discrete approximation for the probability of drawing one of the discrite values, prob_a(:), and one of the corresponding discrete values, a_discr(:)
        !Given the parameters of the ability distribution we derive the support points that gives the best approximation (part B in Kennan)
        DO ii = 1, Nab
           ! cdf values at the fixed support points
           true_cdf_a(ii) = (2d0*DBLE(ii)-1)/(2d0*DBLE(Nab))
           !We compute the inverse of the standard normal at true_cdf_a(ii)
           CALL normal_01_cdf_inv(true_cdf_a(ii), a_discr_temp)
           a_discr(ii) = sigma_a*a_discr_temp + mu_a
        END DO
        
        !We compute the probabilities at the chosen points that best approximate the log-normal distribution (part A in Kennan)
        CALL normal_01_cdf((a_discr(1) - mu_a)/sigma_a,prob_a_temp)
        DO ii = 2, Nab
           CALL normal_01_cdf((a_discr(ii) - mu_a)/sigma_a,prob_a_temp1)
           cdf_a(ii-1) = (prob_a_temp + prob_a_temp1)/2d0
           prob_a_temp = prob_a_temp1
        END DO
        cdf_a(Nab) = 1d0
        
        !We compute the probability of the discrete fine values in our grid
        pability(1) = cdf_a(1)
        DO ii = 2, Nab
           pability(ii) = cdf_a(ii) - cdf_a(ii-1)
        END DO
        
        !We take the exp of a_discr(ii) to make it a log-normal
        a_discr(:) = EXP(a_discr(:))
        
        IF (flag_bootstrap == 0) THEN
           IF(myrank == 0) write(*,'(a)') ''
           IF(myrank == 0) write(*,'(a,3f20.10)') 'true cdf ability', true_cdf_a(:)
           IF(myrank == 0) write(*,'(a,3f20.10)') 'discrete ability', a_discr(:)
           IF(myrank == 0) write(*,'(a,3f20.10)') 'cdf ability', cdf_a(:)
           IF(myrank == 0) write(*,'(a,3f20.10)') 'prob discrete ability', pability(:)
        END IF

        subutil = 0d0
        CALL utilfun(subpcons,qcons_w_riv,theta,uu)
        IF (myrank==0) WRITE(*,*) 'uu,subpcons,subqcons,theta', uu,subpcons,qcons_w_riv,theta
        sumbetadf = 0d0
        DO tt = 1, Tend
           sumbetadf = sumbetadf + betadf**(tt-1)
           subutil(Tend - tt + 1) = uu*sumbetadf
        END DO

        IF (flag_bootstrap == 0 .AND. myrank == 0) THEN
           write(*,*) ''
           write(*,*) 'param', p(:)
           write(*,*) 'electoral probability: alpha_ea, alpha_st, alph_q, alpha_ie', alpha_ea, alpha_st, alpha_q, alpha_ie, alpha_nst
           write(*,*) 'pappeal', pappeal
           write(*,*) 'mu_fine, sigma_fine, delta_fine', mu_fine, sigma_fine, delta_fine
           write(*,*) 'plability', pability
           write(*,*) 'prod fn: pu*ab, pr*ab, educ, const', alphaq(1:4)
           write(*,*) 'prod fn: poly pop', alphaq(5:Nparq)
           write(*,*) 'theta', theta
           write(*,*) 'upower', upower
           write(*,*) 'wages past mayors', alphaw
           write(*,*) 'delta, gamma', delta, gamma
           write(*,*) 'sigma_run, sigma_zzpr', sigma_run, sigma_zzpr
           write(*,*) 'sigma_steal, mu_steal', sigma_steal, mu_steal
           write(*,*) 'cost_run_frac, upower, upower-cost_run_frac', cost_run_frac, upower, upower-cost_run_frac
           write(*,*) 'mu_a, sigma_a', mu_a, sigma_a
           write(*,*) 'beta', beta
           write(*,*) 'eta', eta
           write(*,'(a,11f18.14)') 'param inputs', theta, delta_fine, cost_run_frac, upower, eta, sigma_run, sigma_steal, alpha_ie, alpha_nst, alpha_ea, alpha_q
        END IF
        
        CALL emaxfun()
        
        CALL simulatemun()

        CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)

        !Store the results of each size of perturbation j for each parameter i!
        dsim(i,j,:) = simmom(1:Nmom,1)

        IF (myrank == 0) write(*,'(a,11f20.10)') 'dsim', dsim(i,j,:)

        !----------- COMPUTE CRITERION FUNCTION ------------!

        ! Compute the difference in moments
        diffmom(1:Nmom,1) = simmom(1:Nmom,1) - datamom(1:Nmom,1)

        DO ii=1,Nmom
           temp(ii,1)=SUM(diffmom(:,1)*INV_V(:,ii))
        END DO

        fval=SUM(temp(:,1)*diffmom(:,1))
        IF(myrank ==0) write(*,*) 'i,j,fval',i,j,fval
        fval_vec(j) = fval

        IF (fval < fval_init .AND. i <= np) THEN
           IF (myrank == 0) write(*,'(a,2f20.10)') 'LOWER FVAL', fval, fval_init
           IF (fval < fval_min) THEN
              fval_min = fval
              opt_param = newp
              flag_lowerfval = 1
              IF (myrank == 0) write(*,'(a,2f20.10)') 'NEW LOWER FVAL', fval, fval_min
           END IF
        END IF
        
     END DO
     
     DO j = 2, N_points
        IF(ABS(fval_vec(j) - fval_vec(j-1)) > 0) EXIT
        IF(j == N_points) THEN
           IF(myrank==0) write(*,'(a,1i3,4F20.15)')'No changes in fval',i,fval_vec
        END IF
     END DO
     
     DO j = 1, Nmom

        !Calculating the five point derivative for param i
        dgdp(i,j) = (-25d0*dsim_orig(j)+48d0*dsim(i,1,j)-36d0*dsim(i,2,j)+16d0*dsim(i,3,j)-3d0*dsim(i,4,j))/(12d0*h)

        IF (myrank == 0) write(*,'(a,2i3,7f20.10,1f50.45)') 'np, nmom, dgdp', i, j, dgdp(i,j), -25d0*dsim_orig(j), 48d0*dsim(i,1,j), -36d0*dsim(i,2,j), 16d0*dsim(i,3,j), -3d0*dsim(i,4,j), h, deps
        IF (myrank == 0) write(*,'(a,2i3,6f30.20)') 'np, nmom, dgdp', i, j, dgdp(i,j), -25d0*dsim_orig(j)/(12d0*h), 48d0*dsim(i,1,j)/(12d0*h), -36d0*dsim(i,2,j)/(12d0*h), 16d0*dsim(i,3,j)/(12d0*h), -3d0*dsim(i,4,j)/(12d0*h)
        IF (myrank == 0) write(*,'(a,2i3,1f80.45)') 'np, nmom, den', i, j, 1d0/(12d0*h)

     END DO
     
  END DO

  CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)
  
  DO iparam = 1, np_p
     IF (myrank == 0) write(*,'(a,1i3,11f20.12)') 'Ip', iparam, dgdp(iparam,:)
  END DO

  !Compute the variance-covariance matrix of the estimated parameters
  !First we compute the information matrix Ip
  dg(:,:) = dgdp(1:np,:)
  !Derivative of SMM moments wrt the prod function parameters
  dm(:,:) = dgdp(np+1:np_p,:)

  DO iparam = 1, np
     IF (myrank == 0) write(*,'(a,1i3,23f20.10)') 'Ip', iparam, dg(iparam,:)
  END DO
  DO imom = 1, Nmom
     IF (myrank == 0) write(*,'(a,1i3,23f20.10)') 'inv', imom, INV_V(imom,:)
  END DO

  !VARIANCE COVARIANCE MATRIX IF ONE USES THE EFFICIENT WEIGHTING MATRIX
  !Matrix to be inverted to give the varcov matrix
  DO iparam = 1, np
     DO imom = 1, Nmom
        Vtemp(iparam,imom) = SUM(dg(iparam,:)*INV_V(:,imom))
     END DO
  END DO

  DO i = 1, np
     DO j = 1, np
        VV(i,j) = SUM(Vtemp(i,:)*dg(j,:))
     END DO
     IF (myrank == 0) write(*,'(a,1i3,21f40.20)') 'VV', i, VV(i,:), 1d0/VV(i,i)
  END DO

  !Then we compute its inverse
  Varcov = 0d0
  DO j=1, np
     Varcov(j,j)=1d0
  END DO

  LWORK = np
  IF (ALLOCATED(WORK)) DEALLOCATE(WORK)
  ALLOCATE(WORK(LWORK))
  CALL DPOSV(uplo,np,np,VV,np,Varcov,np,info)

  IF (info /=0 ) THEN
     IF (myrank == 0) write(*,*) 'MATRIX Varcov COULD NOT INVERT', info
  END IF

!!$  Varcov = Varcov/DBLE(Nactmun)
  Varcov = Varcov/SQRT(DBLE(Nactmun))

  DO i = 1, np
     IF (myrank == 0) write(*,'(a,1i3,21f30.15)') 'Varcov', i, Varcov(i,:), 1d0/Varcov(i,i)
  END DO

  DO i = 1, np
     IF (myrank == 0) write(*,'(a,1i3,2f30.15)'), 'i, stderr RATIO', i, pp(i), SQRT(1d0/VV(i,i))/SQRT(DBLE(Nactmun))
  END DO

  DO i = 1, np
     stderr(i) = SQRT(Varcov(i,i))
     IF (myrank == 0) write(*,'(a,1i3,3f25.15)'), 'EFFICIENT: i, param, stderr, tstat', i, pp(i), stderr(i), pp(i)/stderr(i)
  END DO

  ! Part of the variance that accounts for the first-step estimation
  ! We read in the variance matrix of the prod function parameters

  !Then we compute its inverse
  inv_dg = dg

  CALL DGETRF(np, np, inv_dg, np, ipiv, info)

  if (info /= 0) then
      IF (myrank == 0) write(*,*) 'Matrix is numerically singular!', info
  end if

  LWORK = np
  IF (ALLOCATED(WORK)) DEALLOCATE(WORK)
  ALLOCATE(WORK(LWORK))

  CALL DGETRI(np, inv_dg, np, ipiv, work, np, info)

  IF (info /=0 ) THEN
     IF (myrank == 0) write(*,*) 'MATRIX dg COULD NOT INVERT', info
  END IF

  DO i = np+1,np_p
     IF (myrank == 0) write(*,'(a,1i3,20f20.10)') 'dm', i, dgdp(i,:)
  END DO

  DO i = 1, np
     IF (myrank == 0) write(*,'(a,1i3,20f20.10)') 'dg', i, dg(i,:)
  END DO

  DO i = 1, np
     IF (myrank == 0) write(*,'(a,1i3,20f20.10)') 'inv_dg', i, inv_dg(i,:)
  END DO

  ! Check inv_dg is the inverse of dg
  DO i = 1, np
     DO j = 1, np
        Id_temp(i,j) = SUM(dg(i,:)*inv_dg(:,j))
     END DO
     if (myrank == 0) write(*,'(a,11f20.10)') 'Id_temp', Id_temp(i,:)
  END DO
  
  OPEN(25, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/varcov_prod_func_July2023.raw')
  DO j=1, 4
     READ(25,*) temp1(j,:)
     Vpf(j,:) = temp1(j,:)
!!$     IF (myrank == 0) WRITE(*,'(a,12f20.10)') 'Vpf', Vpf(j,:)
     IF (myrank == 0) WRITE(*,'(a,4f20.10)') 'temp1', temp1(j,:)
  END DO
  CLOSE(25) 

  OPEN(25, file='/user/econ/mazzocco/MyPapers/Mayors/Estimation/Data/Data_For_Code/varcov_wage_reg.raw')
  DO j=1, 6
     READ(25,*) temp2(j,:)
     Vw(j,:) = temp2(j,:)
!!$     IF (myrank == 0) WRITE(*,'(a,12f20.10)') 'Vpf', Vpf(j,:)
     IF (myrank == 0) WRITE(*,'(a,6f20.10)') 'temp2', temp2(j,:)
  END DO
  CLOSE(25) 

  Vop = 0d0
  DO j=1, 4
     Vop(j,1:4) = Vpf(j,:)
     IF (myrank == 0) WRITE(*,'(a,1i4,11f16.10)') 'Vop', j, Vop(j,:)
  END DO
  DO j=1, 6
     Vop(4+j,5:10) = Vw(j,:)
     IF (myrank == 0) WRITE(*,'(a,1i4,11f16.10)') 'Vop', j, Vop(4+j,:)
  END DO
    
  Vtemp_op = 0d0
  DO iparam = 1, nop
     DO imom = 1, Nmom
        Vtemp_op(imom,iparam) = SUM(dm(:,imom)*Vop(iparam,:))
     END DO
  END DO

  VV_op = 0d0
  DO i = 1, Nmom
     DO j = 1, Nmom
        VV_op(i,j) = SUM(Vtemp_op(i,:)*dm(:,j))
     END DO
  END DO

  VarCov_temp = 0d0
  DO iparam = 1, np
     DO imom = 1, Nmom
        VarCov_temp(iparam,imom) = SUM(inv_dg(:,iparam)*VV_op(imom,:))
     END DO
  END DO

  VarCov_op = 0d0
  DO i = 1, np
     DO j = 1, np
        VarCov_op(i,j) = SUM(VarCov_temp(i,:)*inv_dg(:,j))
     END DO
  END DO

  VarCov_op = VarCov_op/SQRT(DBLE(Nactmun))
  
  DO i = 1, np
     IF (myrank == 0) write(*,'(a,1i3,11f20.10)') 'Varcov_op', i, Varcov_op(i,:)
  END DO

  TotVarCov = VarCov + VarCov_op

  DO i = 1, np
     IF (myrank == 0) write(*,'(a,1i3,11f20.10)') 'TotVarCov', i, TotVarCov(i,:)
  END DO

  count_statsig = 0
  DO i = 1, np
     IF (myrank == 0) write(*,'(a,1i3,3f25.15)'), 'EFFICIENT TOT: i, param, stderr, tstat', i, pp(i), SQRT(TotVarCov(i,i)), pp(i)/SQRT(TotVarCov(i,i))
     IF (ABS(pp(i)/SQRT(TotVarCov(i,i))) > 1.6d0) count_statsig = count_statsig + 1
  END DO
  IF (myrank == 0) write(*,'(a,1i5)') 'NUMBER OF SIGNIFICANT PARAMETERS', count_statsig
!!$  IF (count_statsig >= np) THEN
!!$     IF (myrank == 0) write(*,*) 'ALL PARAMETERS SIGNIFICANT'
!!$  END IF
  
  VarCovPol = VarCov
  DO i = 1, np
     VarCovPol(i,i) = TotVarCov(i,i)
  END DO
  
  IF (myrank == 0) THEN

     ! We save the estimated parameters and their variance matrix
     OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Data/Data_For_Code/estparam_stderrs_tot.dat", STATUS='REPLACE')
     WRITE(28,'(11f25.18)') pp(1:np)
     WRITE(*,'(11f25.18)') pp(1:np)
     DO i = 1, np
        WRITE(28,'(11f25.18)') TotVarCov(i,:)
        WRITE(*,'(11f25.18)') TotVarCov(i,:)
     END DO
     WRITE(*,*) ''
     CLOSE(28)
     
     ! We save the estimated parameters and the variance matrix with the standard errors in the main diagonal
     OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Data/Data_For_Code/estparam_stderrs_pol.dat", STATUS='REPLACE')
     WRITE(28,'(11f25.18)') pp(1:np)
     DO i = 1, np
        WRITE(28,'(11f25.18)') VarCovPol(i,:)
        WRITE(*,'(11f25.18)') VarCovPol(i,:)
     END DO
     WRITE(*,*) ''
     CLOSE(28)

     ! We save the estimated parameters and the variance matrix without considering the first step estimates
     OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Data/Data_For_Code/estparam_stderrs.dat", STATUS='REPLACE')
     WRITE(28,'(11f25.18)') pp(1:np)
     DO i = 1, np
        WRITE(28,'(11f25.18)') VarCov(i,:)
        WRITE(*,'(11f25.18)') VarCov(i,:)
     END DO
     CLOSE(28)
     
  END IF
  
END SUBROUTINE standarderrors
